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ABSTRACT 

Motivated by the recent detection of a large number of embedded young stellar objects (YSOs) with 
mass acc retion rat es that are inconsistent with the predictions of the standard model of inside-out 
collapse ([ShulllQTTi) . we perform a series on numerical hydrodynamic simulations of the gravitational 
collapse of molecular cloud cores with various initial masses, rotation rates, and sizes. We focus on 
the early Class I stage of stellar evolution when circumstellar disks are exposed to high rates of mass 
deposition from infalling envelopes. Our numerical modeling reproduces the large observed spread in 
accretion rates inferred for embedded YSOs in Perseus, Serpens, and Ophiuchus star forming regions by 
lEnoch et al] (|2009f ). yielding 37%-75% of objects with "sub-Shu" accretion rates M < 10"^ M© yr"! 
and l%-2% of objects with "super-Shu" accretion rates M > 10~^ Mq yr~^. Mass accretion rates 
in the Class I stage have a log-normal distribution, with its shape controlled by disk viscosity and 
disk temperature. The spread in M is greater in models with lower viscosity and smaller in models 
with higher viscosity and higher disk temperature, suggesting that gravitational instability may be a 
dominant cause of the observed diversity in M in embedded YSOs. Our modeling predicts a weak 
dependence between the mean mass accretion rates and stellar masses in the Class I stage, in sharp 
contrast to the corresponding steep dependence for evolved T Tauri stars and brown dwarfs. 

Subject headings: circumstellar matter — planetary systems: protoplanetary disks — hydrodynamics 
— ISM: clouds — stars: formation 



1. INTRODUCTION 

It has recently become evident that young stellar ob- 
jects (YSOs) in the embedded phase of stellar evolution 
exhibit a variety of mass accretion rates, which range be- 
tween M ^ 10"^ M© yr~^ and M 10" '^ Mq yr'^ (e.g. 
iDunham etld] l2006HEnoch et all 1200 9^ and ofte n can- 
not be ac counted fo r in the standard model of inside-out 
collapse (jShulllOTTl ). According to this model, the infall 
rate of gas from an envelope onto a star/disk system is 

Miniedi ~ c^/G, where Cg is the sound speed. For the typi- 
cal gas temperatures in star forming regions rg=10-20 K, 
the standard model predicts infall rates of order 1.2 - 
4.2x10-6 Mq yr-i. According to lEnoch et aH (|2009f ). 
however, more than 50% of YSOs in Perseus, Serpens, 
and Ophiuchus star-forming regions have inferred mass 
accretion rates 10"^ Mq yr"^ < M < 10"^ Mq yr'^ 
which are considerably lower than the infall rates. At 
the same time, the envelope m ass decreases abou t as 
predicted by the standard model ([Enoch et al.ll2009f ) , im- 
plying that the observationally inferred and theoretically 
derived infall rates Minfaii agree with each other. 

The obvi ous mismatch between M and Minfaii was first 
noticed by iKenvon et al.l (|1990l ) based on the compar- 
ison between typical lifetimes of embedded YSOs and 
inferred bolometric luminosities — the accretion luminos- 
ity ~ 10~6 — 10~^ Mq yr~^ inferred from the dura- 
tion of the embedded phase a few x 10^ yr turned out 
to be much larger than the typical bolometric luminos- 
ity ~ 10~^ Mq yr~^. This "luminosity problem" can 
be solved by the infalling material first piling up in a 
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circumstellar disk and then being accreted onto a pro- 
tostar mostly in short-lived FU-Ori-like episodes with 
accretion rates M > 10^^ Mq yr~^. Numerous mech- 
anisms for this episodic accretion have been proposed 
in the pa st. These include the t hermal ionization in- 
stabi li tv (iLin fc Papaloizoul 119851 : iHartmann fc KenvonI 
IT985I: iBell fc LinI Il994fl. cl ose e ncounters in binary 
systems ( Bonnell fc Basteinl Il992l ). gravitational insta- 
bility and fragmentation in embedded circumstellar 
disks (jVorobvov fc B as"?! 2OO51 I2OO60 . a combination the 
magneto-rotational instability in the inner disk regions 
and gravitationa l instability in the outer disk regions 
(IZhu et al.ll2009f). and close encounters in young stellar 
clusters (jPfalzner et al.ll2003^ . 

Most theoretical and numerical studies of episodic ac- 
cretion have been focused so far on explaining observ- 
able characteristics of the FU-Ori-like accretion bursts. 
In this paper, however, we mostly focus on the quiescent 
phase of accretion between the bursts. This phase has re- 
cently gained much interest due to a discovery of very low 
luminosity objects or VELLOs in dense molecular cores 
that have previou s ly been thought of as "starless" (e.g. 
Young et al]|2002t iKauffmann et afl 120051 : [Bolarke et al.l 



200a IStecklum et al.ll2007t ). One possible explanation 



for their low luminosity {L < 0.1 Lq) is a quiescent 
accretion, in which M has to be more than an or- 
der of magnitude lower than the Shu infall rates 1.2- 
4.2xl0~6 Mq yr~^, suggesting again that the protostel- 
lar accretion history may not be as simple as predicted 
by the standard model. 

In this paper, we provide a comprehensive study of 
the mass accretion rates M and disk infall rates Minfaii 
in the embedded phase of stellar evolution (focusing 
mostly on the Class I stage). We use numerical hydro- 
dynamic simulations of the cloud core collapse with a 
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self-consistent disk formation and long-term evolution. 
We run a large set of model covering a wide range of 
cloud core masses and angular velocities. The general 
problem of explaining the large spread of observed pro- 
tostcUar accretion rates in the Class I stage is addressed 
using ideas from the burst mode of accret ion discovered 
recently bv lVorobvov k Basul (|2005l . l2006f) . 

2. DESCRIPTION OF NUMERICAL MODEL 

Our numerical mod e l is sim ilar to that used recently 
by IVorobvov &: Basul ()2009aD to simulate the secular 
evolution of viscous and self-gravitating circumstellar 
disks and the mass accretion rates in the late evolution 
phase of T Tauri stars (T TSs) and brown dwarfs (BDs) 
(jVorobvov fc Basul[200&bl ). For the reader's convenience, 
we briefly review the basic concept and equations. 

We employ numerical hydrodynamic simulations in the 
thin-disk approximation to compute the evolution of ro- 
tating, gravitationally unstable cloud cores with various 
initial masses Md and ratios of rotational to gravita- 
tional energy (3. We start our numerical integration in 
the pre-stellar phase, which is characterized by a col- 
lapsing starless cloud core, continue into the embedded 
phase of stellar evolution, which sees the formation of 
a star/disk/envelope system, and terminate our simula- 
tions in the late accretion phase, when most of the cloud 
core has accreted onto the star/disk system. Once the 
disk has self-consistently formed, it occupies the inner- 
most regions of our numerical grid, while the infalling 
cloud core (the envelope) occupies the rest of the grid. 
This ensures that the mass infall rate onto the disk Minfaii 
is accurately determined by the dynamics of the gas in 
the envelope. 

The basic equations of mass and momentum transport 
in the thin-disk approximation are 
r)y 

— = -V,.(S^,), (1) 

fin 

S^ = "Vp7' + Eg, + (V.n),, (2) 

where E is the mass surface density, V = Pdz is 
the vertically integrated form of the gas pressure P, 
Z is the radially and azimuthally varying vertical scale 
height, Vp = VrV + v^(f) is the velocity in the disk plane, 
9p — drf" + g(t>4' is the gravitational acceleration in the 

disk plane, and Vp = rd/dr + (f)r~^d/d(j) is the gradient 
along the planar coordinates of the disk. The gravita- 
tional acceleration includes both the gravity of a cen- 
tral point object (when formed) and the self-gravity of a 
circumstellar disk and envelope. The latter component 
is found by solving the Poisson integral using the convo- 
lution theorem. The viscous stress tensor 11 is expressed 
as 

n = 2Ei/^Vw-i(V-w)e^ , (3) 

where Vv is a symmetrized velocity gradient tensor, e 
is the unit tensor, and v is the kinematic viscosity. The 
compon ents of (V • II)c in polar coo rdinates (r, cj)) can be 
found in IVorobvov fc 13asul (|2009al) . 

The best candidate for viscosity in circumstellar disks 
is turbulence induced by the magneto-rotational insta- 
bility (MRI), though other mechanisms such as non- 
linear hydrodynamic turbulence cannot be completely 



eliminated due to the large Reynolds numbers involved. 
We make no assumptions about the source of viscos- 
ity and para meterize its magnitude us ing the usual a- 
prescription (jShakura fc SunvaevlllQTl ) 

v = ac.sZ, (4) 

where = dV/dYj is the effective sound speed of (gen- 
erally) non- isothermal gas. The vertical scale height Z is 
determined in every computational cell and at every time 
step of integration using an assumption of local hydro- 
static equilibrium in t he gravitational field of th e central 
star and the disk (see IVorobvov fc Basull2009ah . 

For most of the numerical simulations in this paper, we 
us e a = 0.01. This ch o ice is m otivated by the recent work 
of IVorobvov fc Baiul (|2009al) . who studied numerically 
the secular evolution of viscous and self-gravitating disks. 
They found that if circumstellar disks around solar-mass 
protostars could generate and sustain turbulence than 
the temporally and spatially averaged a should lie in the 
10^3 _ iQ-2 lijnits. Smaller values of a (< lO"'^) have 
little effect on the resultant mass accretion history, while 
larger values {a > 10~^) destroy circumstellar disks dur- 
ing less than 1.0 Myr of evolution and are thus unlikely 
from the point of view of disk longevity. We note that we 
have intentionally taken the largest possible value for a, 
since we want to assess the maximum effect that viscos- 
ity may have on the accretion history in the embedded 
phase. The case of a smaller a = 10~^ will be consid- 
ered in brief in Section 15721 We also note that viscosity is 
introduced in numerical simulations only after disk for- 
mation. In the pre-disk phase, the a-parameter is set to 
zero. 

Equations ([T]) and ([2]) are closed with a barotropic 
equation that makes a smooth transition from isother- 
mal to adiabatic evolution at E = = 36.2 g cm^^: 

V^cl^ + cl^^l^^y , (5) 

where Cs = 0.188 km s~^ is the sound speed in the begin- 
ning of numerical simulations (corresponding to the ini- 
tial gas temperature of 10 K) and 7 = 1.4. The adopted 
value for Ecr corresponds to the gas volume number den- 
sity of 10~^^ cm~'^ for a disk in the vertical hydrostatic 
equilibrium at temperature 10 K. The effect of larger 7 
and, as a consequence, hotter circumstellar disks is ex- 
plored in Section 15.31 

Equations ([T]) and ([2]) are solved in polar coordinates 
(r, (/)) on a numerical grid with 128 x 128 points. We have 
found that an increase in the resolution to 256 x 256 
grid zones makes little influence on the accretion his- 
tory but helps to save a considerable amount of CPU 
time and consider many more model cloud cores. We 
use the method of finite differences with a time-explicit, 
operator-split solution procedure similar in methodology 
to the ZEUS code. Advection is performed using the 
second-order van Leer scheme. The radial points are log- 
arithmically spaced. The innermost grid point is located 
at — 5 AU, and the size of the first adjacent cell varies 
in the 0.17-0.36 AU range depending on the cloud core 
size. We introduce a "sink cell" at r < 5 AU, which 
represents the central star plus some circumstellar disk 
material, and impose a free inflow inner boundary condi- 
tion. The outer boundary is reflecting. A small amount 
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of artificial viscosity is added to the code, though the as- 
sociated artificial viscosity torques are always negligible 
in comparison with the gravitational torques. 

3. INITIAL CONDITIONS 

We start our numerical simulations from starless cloud 
cores, which have surface densities S and angular veloci- 
ties n typical for a collapsing axisymmetric magnetically 
supercritical core (|Basulll997D 




where fJo is the central angular velocity, tq is the ra- 
dius of central near-constant-density plateau defined as 
To = fcc^/(GI]o) and k = These initial profiles 

are characterized by the important dimensionless free 
parameter j] = ilprg/c^ and have the property that the 
asymptotic (r ^ rg) ratio of centrifugal to gravitational 
acceleration has magnitude v^??- The centrifugal radius 
of a mass shell initially located at radius r is estimated to 
be Tcf = p /{Gni) — V2r]r, where j = ilr^ is the specific 
angular momentum. We note that rj is similar in mag- 
nitude to the ratio of rotational to gravitational energy 
P = -Erot /-E'grav , whcrc the rotational and gravitational 
energies are calculated as 

rout rout 

E-cot — 2i: j racY^rdr, i?g,.av = — Stt j rg^^rdr. 

. (8) 

Here, Qc = r is the centrifugal acceleration, and rout 
is the outer cloud core radius. The numerical relation- 
ship between the two parameters is /3 — O.9I77. The gas 
has a mean molecular mass 2.33 mn and cloud cores are 
initially isothermal with temperature T = 10 K. 

We present results from nine sets of models, the pa- 
rameters of which are detailed in Table [T] Each set is 
characterized by a distinct ratio f3 of rotational to grav- 
itational energy. The adopted values of f3 lie within the 
limits inferred bv ICaselli et al.l ()2002D for dense molecu- 
lar cloud cores, /3 = (10"'' — 0.07). The product ro^o for 
every model in a given set is kept constant to enforce the 
P = const condition (the initial sound speed Cg — 0.188 
km s~^ is equal for all models). In addition, the ratio 
?'out/''o is set to 6.0 to generate gravitationally unstable 
truncated cores of similar form. 

The actual procedure for generating a specific cloud 
core with a given value of p is as follows. First, we 
choose the outer cloud core radius rout and find rp from 
the condition rout/j'o = 6-0. Then, we find the central 
surface density Eq from the relation tq = ^/2c^/(ttGT,q) 
and determine the resulting cloud core mass Afci from 
Equation ([6]). Finally, the central angular velocity fto is 
found from the condition /3 = O.9I77 = 0.91f2§rQ/Cs. 

In total, we have simulated numerically the time evo- 
lution of 78 cloud cores. We note that model sets with 
greater /3 have on general a smaller number of mod- 
els due to numerical difficulties associated with mod- 
eling of massive cloud cores with high angular veloc- 
ities. The resulting initial cloud core mass function 



dN/dAIci = M™ has exponents mi = -0.4 ± 0.1 in the 
0.1 Mq < Mci < 1.0 Mq mass range and m2 = -1.3±0.2 
in the 1.0 Mq < Mc\ < 3.0 Mq range. These values are 
somewhat shallower than usually inferred f rom nearby 
star-forming regions (e.g. lEnoch et al.l |2006'). The slope 
of the cloud core mass function varies in different star 
forming regions and the effect of this variation will be 
considered in a follow-up study. 

4. TIME EVOLUTION OF MASS ACCRETION RATES 

In this section, we review the accretion history in the 
early phase of stellar evolution as derived using our nu- 
merical hydrodynamic simulations. We calculate the in- 
stantaneous mass accretion rate M = — 27rrsci'rao^(^sc) 
as the mass passing through the sink cell per one time 
step of numerical integration (which in physical units is 
usually equal to 10-20 days). We note that the size of the 
sink cell rsc = 5 AU is larger than the stellar radius. The 
inner disk at r < 5 AU may add additional variability to 
the accretion rates, in particular due to the th ermal ion- 
ization instability (jBell fc Lid 1 1994( 1 or MRI (jZhu et all 
I2009f ). These effects may somewhat alter the temporal 
behaviour of the actual accretion rates onto the stellar 
surface. Protostellar jets may also decrease the rate of 
mass deposition onto the star by about 10%. We leave 
these complicated phenomena for a future study. 

Figure [1] shows the time behaviour of M (solid lines) 
for 12 models characterized by distinct initial cloud core 
masses Mc\ (in Mq) and ratios of the rotational to grav- 
itational energy /?. In addition, the dashed lines present 
mass infall rates calculated as Minfaii = — 27rrWrS(r), 
where r = 600 AU. The disk radius in our models rarely 
exceeds 600 AU and Minfaii represents the rate of gas 
deposition from an infalling envelope onto a disk. The 
vertical dotted lines mark the onset of the Class I (left) 
and Class II (right) stages of stellar evolution, as inferred 
from the mass remaining in the e nvelope. Following a 
prescription of lAndre et al.l ([Till) (which is most use- 
ful for numerical simulations), we assume that the Class 
I/Class II stages ensue when 50%/90% of the initial cloud 
core mass, respectively, has been accreted onto the star 
plus disk system. Other classification schemes may in- 
troduce a systematic shift of order unity. 

The early evolution of the mass accretion rate is sim- 
ilar in all models — M reaches a maximum value of ^ 
2 X 10~^ Mq yr~^ and then settles at a near-constant 
value of order of 10~^ Mq yr~^. This early phase of 
near-constant accretion ends when the first layer of in- 
falling material hits the centrifugal radius at 5.0 AU and 
a disk starts to form in our numerical simulations. The 
subsequent evolution of mass accretion is controlled by 
the physical processes in the disk and the rate of mass 
loading Afinfaii onto the disk. A continuing deposition 
of matter from the infalling envelope, as indicated by 
the dashed lines in Figure [1] leads to the development of 
gravitational instability and fragmentation in the disk. 
Fragments typically form at r > 40 AU, have typical 
masses of up to 10-20 Jupiters, sizes of several AU, num- 
ber densities of up to 10^^ — lO^'* cm~^, and are pressure 
supported against their own gravity. These fragments are 
quickly driven into the inner disk regions via the gravi- 
tational exchange of angular momentum with the spiral 
arms and trigger a short-lived burst of mass accretion 
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TABLE 1 

Model parameters 



Note. 



Set 


/3 














N 


1 


1.18 X 10"^ 


0.27 - 1.00 


1380 


- 5180 


(0.8 


- 3.0) X 10'' 


0.8 - 3.0 


8 


2 


2.29 X 10-^ 


0.37 - 2.24 


860 


- 5180 


(0.5 


- 3.0) X lO'' 


0.5 - 3.0 


11 


3 


3.37 X 10-3 


0.45 - 3.40 


690 


- 5180 


(0.4 


- 3.0) X 10* 


0.4-3.0 


14 


4 


5.10 X 10-3 


0.53 - 6.70 


410 


- 5180 


(0.2 


- 3.0) X 10'' 


0.24 - 3.0 


12 


5 


7.25 X 10-3 


0.67 - 8.30 


414 


- 5180 


(0.2 


- 3.0) X 10" 


0.24 - 3.0 


9 


6 


1.20 X 10-2 


1.20 - 17.0 


240 


- 3450 


(0.14 


- 2.0) X 10" 


0.14 - 2.0 


11 


7 


1.42 X 10-2 


2.30 - 14.0 


340 


- 2070 


(0.2 


- 1.2) X 10" 


0.2 - 1.2 


4 


8 


2.00 X 10-2 


1.60 - 22.9 


240 


- 3450 


(0.14 


- 2.0) X 10" 


0.14 - 2.0 


5 


9 


2.90 X 10-2 


2.00 - 20.0 


340 


- 3450 


(0.2 


- 2.0) X 10" 


0.2 - 2.0 


4 


5 are 


in AU, angular 


velocities in 


km s ^ pc ^ , 


masses 


in Mq, and A'' 


is the number of 
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Fig. 1. — Mass accretion rates M (solid lines) and mass infall rates A/infaii (dashed lines) versus time elapsed since the beginning of 
numerical simulations. Twelve models with distinct initial cloud core masses M^i (in Mq) and ratios of rotational to gravitational energy 
f3 are presented. In particular, the top/middle/bottom rows have /3 = 2.0 X 10—2/1.2 X 10-2/5.1 X IQ- 3, respectively. The vertical dotted 
lines mark the onset of the Class I stage (left) and Class II stage (right) of stellar evolution. 



when passing through the sink cell. This phenomenon is 
called th e burst mode of accretion and is investigated in 
detail bv l\^robvov fc Basul (|2005l . [20061) . The ultimate 
fate of the fragments is uncertain and is largely depen- 
dent on how quickly they can contract from their initial 
size of several AU to the planetary size to a void tidal 
destruction. According to iHelled et al.l (|2006f ). the con- 
traction time for a Jupiter-mass clump to reach a central 
temperature of 2000 K, i.e., the temperature required to 
dissociate H2 to trigger rapid collapse, is 3 x 10^ yr. How- 



ever, considering a fast time scale of inward radial migra- 
tion to the inner 5 AU — a few tens of orbital periods at 
most"^ — we believe that most of these fragments will be 
tidally destroyed when approaching the central star, thus 
converting its gravitational energy to the accretion lumi- 
nosity and producing an FU-Ori-like luminosity burst. 
Figure [1] demonstrates that the burst phenomenon is 

3 The animation of this process can be viewed at 
www.astro.uwo.ca/ ~vorobyov 
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Mass accretion rate (M@ yr"^) 

Fig. 2. — The normalized number of Class_I YSOs versus mass 
accretion rate as inferred by Enoch et al. ( 2009!) ™ Perseus, Serpens 
and Ophiuchus star forming regions. The solid and dashed lines 
are the log-normal and gaussian fits to the observed data. 

more pronounced in models with greater initial cloud 
core masses Md- Indeed, more massive cloud cores 
have larger outer radii rout and, as a consequence, larger 
centrifugal radii rd- It then follows that models with 
greater Md are expected to form more massive disks, 
which are easier susceptible to gravitational instability 
and fragmentation. Another important feature that can 
be seen in Figure [1] is that the bursts cease when the 
rate of mass deposition onto the disk Afinfaii (dashed 
lines) falls below that of the mass accretion M (solid 
lines). This clearly demonstrates the destabilizing influ- 
ence that the infalling envelope has on the disk. The 
burst phenomenon is mostly localised to the Class and 
Class I stages of stellar evolution, with only a few bursts 
taking place in the early Class II stage. 

5. ACCRETION VARIATIONS IN THE EMBEDDED PHASE 
OF STELLAR EVOLUTION 

Accretion variations appear to be an intrinsic property 
of self-gravitating disks in the early embedded phase of 
stellar evolution. Young circumstellar disks, when ex- 
posed to mass deposition from infalling envelopes, un- 
dergo repeating cycles of gravitational destabilization 
and fragmentation. Prolonged periods of rather low 
accretion M < 10~^ Mq yr~^, when the disk may 
be gravitationally unstable but is yet stable to frag- 
mentation, are interspersed with short-lived bursts with 
M > 10~^ Mq yr~^, when the disk forms fragments 
and drives them onto the central star. It is therefore 
interesting to directly compare our predicted spread in 
M with that inferred in YSOs in P erseus, Serpens, an d 
Ophiuchus star-forming regions bv lEnoch etall (|2009D . 
In the remaining text, we focus on the Class I stage of 
stellar evolution, leaving the Class stage for a follow-up 
study. 

5.1. Distribution Functions of Accretion Rates 

We start by constructing the distribution function 
(DF) of accretion rates, which shows the normalized 
number o f YSOs versus accr etion rate using data ob- 
tained by lEnoch et al.l (|2009f ) in Perseus, Serpens, and 
Ophiuchus star-forming regions. We distribute the ob- 



served accretion rates in 20 logaritmically spaced bins 
between 10"^ Mq yr~i and 10"'' Mq yr~^ The normal- 
ization is done by first finding a bin that happens to have 
the maximum number of YSOs and then normalizing all 
bins to this maximum number. The resultant normal- 
ized distribution is shown in Figure O The dashed/solid 
lines show gaussian/log-normal fits to the observed data. 
Neither gaussian nor log-normal functions yield good re- 
sults, though the latter function is more preferable than 
the former. There is a lack of well defined peak in the 
observed accretion rates, which may to some extent be 
attributed to a small number of observed objects. The 
log-normal fit yields the following relation between the 
normalized number of YSOs Nygo and observationally 
inferred accretion rates M 



YSO 



0.7 exp 



ln(M/Mo) 
2^1 



(9) 



where Mq = 4 x 10^^ Mq y r~^ is the ge o metri c mean 
accretion rate. According to lEnoch et al.l (|2009D , there 
are approximately 55% of YSOs with sub-Shu accretion 
rates {M < 10"^ yr^^) and about 5% of YSOs with 
super-Shu rates (M > 10"^ Mq yr~^). 

One feature in Figure |5| is worth special attention — 
there is a sizeable portion of YSOs with accretion 
rates smaller than 10"'' Mq yr~^. These objects 
are usually characterized by bolometric luminosity 
Lhoi < 0.1 Lq and constitute a group of objects 
called VELLOs or very low luminosity objects (e.g. 
Young et al."2 002l: iKaufaiann et all 120051 : iBourke et al.l 
2006; Stecklum et al.ll2007n. Oii t of 89 Class I YSOs in 
the compilation of lEnoch et all (|2009[ ). 19 objects are 
VELLOs, which is roughly 21%. 

To compare our models with observations, we calculate 
M every 20 yr in the Class I stage of stellar evolution. 
In the following text, we refer to these calculations as 
" accretion measurements" or simply " measurements" by 
analogy to observations. We note that due to the use 
of the sink cell in our numerical simulations, the Class I 
stage in some models may start before the disk begins to 
form (see Fig. [T]). In such rare cases, accretion rates are 
calculated only starting from the disk formation epoch. 
Next, we construct the distribution function of accretion 
rates showing the normalized number of accretion mea- 
surements versus the accretion rate in 50 logarithmically 
spaced bins between 10~* Mq yr~^ and lO^'* Mq yr"^. 
The normalization is done by first finding a bin that hap- 
pens to have the maximum number of accretion measure- 
ments and then normalizing all bins to this maximum 
number. 

Figure |3| presents DFs in the Class I stage for the 
same set of models as in Figure [TJ Model accretion 
rates approximately follow a log-normal distribution, the 
geometric mean values are mostly located in the 0.5- 
1.0x10"^ Mq yr~^ range. There is a mild tendency for 
models with larger values of (3 to have a wider range of 
accretion rates. It is also seen that the magnitude of 
intrinsic variability in most models is insufficient to ac- 
count for a wide spread in observed mass accretion rates 
(Fig. [2]), implying that some object-to-object variations 
(in the initial cloud core masses and/or rotation rates) 
are needed to reproduce the observed spread. 
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Fig. 3. — Distribution functions of accretion rates showing the normalized number of accretion measurements in the Class I stage for 12 
models with the same parameters as in Figure [T] 



Another interesting feature in Figure [3] is a clear trend 
that, for the same (3, the peak M decreases with in- 
creasing cloud core mass. This tendency is especially 
apparent for cloud cores with sufficiently high rotation 
rates, /3 > 1.2 x fO"^ ^nd can be explained by a growing 
strength of the burst mode of accretion along the se- 
quence of increasing cloud core masses. Indeed, Figure [T] 
demonstrates that, for the same /3, the number of ac- 
cretion bursts increases with increasing cloud core mass, 
which is explained by growing disk masses and increasing 
propensity of disks to fragmentation. At the same time, 
periods of quiescent accretion between the bursts become 
also more frequent and are characterized by lower rates. 
This acts to decrease the peak M along the sequence of 
increasing cloud core masses. By the same reason, the 
peak M has a weak dependence on f3 for cloud cores of 
similar mass. 

The fact that M in our numerical simulations are log- 
normally distributed is interesting. It suggests that the 
accretion process is not statistically independent but it 
has a memory effect — short periods of elevated accretion 
are, as a rule, followed by prolonged periods of quiescent 
accretion and vice versa. This pattern of time behaviour 
is not chaotic but is rather controlled by disk physics, 
implying a causal link between accretion events that are 
separated not too far in time. 



As the next step, we want to construct the integrated 
DF of accretion rates, which would include data from 
all models in Table [T] Every individual model is char- 
acterized by a distinct duration of the Class I stage, 
which increases along the sequence of increasing stellar 
masses. This means that models with longer lifetimes 
of the Class I stage are more significant statistically and 
have a larger number of accretion measurements than 
models with shorter Class I lifetimes. At the same time, 
objects at the lower mass end (with shorter Class I life- 
times) are expected to be more abundant. The letter 
tendency, however, may be counterbalanced by accretion 
measurements biased toward more luminous (and mas- 
sive) objects. Therefore, we first sum up non-normalized 
DFs for every models and then normalize the integrated 
DF in the same manner as described above, thus accen- 
tuating models with a larger number of accretion mea- 
surements. 

The resultant integrated DF is presented in the upper- 
left panel of Figure [H The solid line shows a log-normal 
fit to the model data, while the dashed line is a gaus- 
sian fit. It is evident that the log- normal function yields 
a much better fit to the model data than the gaussian 
one. In particular, we find the following relation be- 
tween the normalized number of accretion measurements 
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-^modoi and the mass accretion rates M 



N„ 



Ddol 



0.97 exp 



ln(Af/Afo) 
0J5 



(10) 



where Mq = 8-8 x 10~^ Mq yr~^ is the geometric mean 
value. Our modehng predicts that about 45 ± 5% of 
YSOs are expected to have "sub-Shu" accretion rates 
with M < 10~^ Mq yr~^. Approximately the same 
number of YSOs are expected to have accretion rates 
that are roughly consistent with the standard model 
10^6 Mq yr-i < M < 10'^ Mq yr'^, and only 2% 
of YSOs have "super-Shu" accretion rates with M > 
10^^ Mq yr~^. These numbers are in far e agreement 
with what obtained by lEnoch et al.l (|2009( l. However, 
our numerical modeling predicts a much smaller frac- 
tion of VELLOs with ^ M < 10"^ M q yr'^, - 1.0%, in 
contrast to ~ 



21% in lEnoch et al.l ([2009). This is most 
likely related to the fact we have adopted too large a 
value for the a-parameter, a = 0.01, thus overempha- 
sizing the effect of viscosity in the early disk evolution. 
It is known that viscosity acts to suppress gravitational 
instability and associated large var iations in the accre- 
tion rates (jVorobvov fc Basull2009al ). reducing the num- 
ber and frequency of FU Orionis or VELLO events. The 
effect of a smaller a-parameter on our model mass accre- 
tion rates is discussed below. 

5.2. The effect of lower disk viscosity 

We have considered so far model disks with a = 10~^. 
As discussed in Section [2l this value of the spatially 
and temporally averaged a-parameter is a realistic up- 
per limit, though some short-living fluctuations towards 
larger values are possible. In this section, we consider 
disks with an order of magnitude lower a-parameter. The 
upper-right panel in Figure|4|presents the integrated dis- 
tribution function of mass accretion rates for all models 
from Table [1] with a set to 10"'^. The solid line is the 
log-normal fit to the model data described as 



iVmodci = 0.94 exp ( -i 



ln(Af/Mo) 
L25 



(11) 



where Mq = 3.8 x lO"'^ Mq yr^^. 

Numerical simulations with a = 10^'^ reveal a sub- 
stantial number (~ 15%) of VELLOs with M < 
10~^ Mq yr~^, in f are ag reement with Figure ID and 
data of lEnoch et al.l (|2009f) . On the other hand, mod- 
els with a = 10-3 predict that about 74 ± 4% of YSOs 
have sub-Shu accretion rates with M < 10~^ Mq yr"^ 
and only 24 ± 4% of them have accretion rates that are 
consistent with the standard model 10"^ Mq yr"^ < 
M < 10-^ Mq yr-^ in contrast to 55% and 40%, re- 
spectively, found by lEnoch et al.l (120091 ). It is evident 
that the a = 10^^ models somewhat overestimate the 
number of sub-Shu objects as inferred from observations. 
This suggests that a preferable value for the temporally 
and spatially averaged a-parameter may lie somewhere 
between 10"'^ and 10~^. In both limiting cases, about 
2% of YSOs are expected to have super-Shu accretion 
rates with M > 10~^ Mq yr~^ 



5.3. The effect of higher disk temperature 

In this section, we briefly consider the effect that a 
greater ratio of specific heats 7 = 1.6 may have on the 
mass accretion rates. An increase in 7 in barotropic mod- 
els leads to an increase in disk temperature. For instance, 
in the 7 = 1.4 case a mean disk temperature at r=10 AU 
(in the Class I stage) for nine models of Figure [1] is 29- 
32 K, while in the 7 = 1.6 case it equals to 45-50 K. The 
corresponding integrated distribution function of accre- 
tion rates is shown in the lower-left panel of Figure |4l 
The a-parameter is set to 10~^. As usual, the solid line 
is the log-normal fit to the model data described by the 
following relation 



model 



1.0 exp 



1 



ln(M/Mo) 
068 



(12) 



where Mq = 1.4 x 10 ® Mq yr ^. A visual com- 
parison of the upper-left and lower-left panels demon- 
strates that the spread in M in the 7 = 1.6 models 
is narrower than in the 7 = 1.4 case, which is a di- 
rect consequence of weaker gravitational instability in 
warmer disks. More specifically, there are 37% ± 5% 
sub-Shu objects and 62% ± 5% "standard" objects with 
10"^ Mq yr-i < M < 10"^ Mq yr'^, whereas ob- 
servational data indicate a prevalence of sub-Shu accre- 
tors (~ 55%). The shortage of sub-Shu accretors in the 
7 = 1.6 case as compared to Eno ch et al.l |2009) suggests 
that the observed circumstellar disks may be relatively 
cold. Alternatively, a combination of warmer disk tem- 
perature and lower turbulent v iscosity (with a_~ lO^'^) 
may also provide a better fit to lEnoch et al.l (|2009[ ). 

5.4. Purely viscous disks 

As the last numerical experiment, we consider cir- 
cumstellar disks in which mass and angular momentum 
transport is governed exclusively by turbulent viscos- 
ity. This is done with the purpose to investigate the 
role of disk self-gravity in the early embedded phase of 
stellar evolution. We artificially set disk self-gravity to 
zero immediately after disk formation, but the gravity 
of the central star is kept unchanged. The resultant in- 
tegrated DF of accretion rates for all models of Table [T] 
is shown in the lower-right panel of Figure [H Purely 
viscous models greatly underestimate the number of ob- 
jects with "standard" accretion rates 10~^ Mq yr^^ < 
M < 10"^ Mq yr~^. In particular, modeling predicts 
7 % of such objects, i n contrast to ~ 45% in the sample 
of lEnoch et all (|2009f ). There are no objects with super- 
Shu accretion rates M > 10~^ Mq yr~^. This example 
serves to demonstrate the importance of disk self-gravity 
in the early disk evolution. 

6. ACCRETION RATES VERSUS STELLAR MASSES 

It has recently been noticed that the mass accretion 
rates of TTSs and BDs of 0.5-3.0 Myr age have a strong 
dependence on the central object mass Af* of the form 
M PC with n « 2.0 (see e.g. iMuzerolle erahl 120031: 
iNatta et al.ll2003 l. A numerical model that is used in 
the p resent work can in princ ipal reproduce this rela- 
tion (jVorobvov fc Basull2009b[ ). predicting a somewhat 
steeper dependence for BDs and low-mass TTSs {n w 
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Fig. 4. — Integrated distribution function of accretion rates in the Class I stage showing the normahzed number of accretion measurements 
^modcl versus mass accretion rate M for all models of Table [T] Every panel is characterized by distinct values of the o-parameter and the 
ratio of specific heats 7 as indicated. In addition, the lower-right panel has disk self-gravity artificially set to zero. The solid and dashed 
lines are the log-normal and gaussian fits to the model data, respectively. 

2.9 ± 0.5), as observed, and a shallower dependence for 
intermediate- and upper-mass TTSs {n « 1.5 ± 0.1), as 
also observed. It is therefore very interesting to see if a 
similar relation is expected in the early embedded phase 
of stellar evolution. 

For every model from Table [TJ we calculate the ge- 
ometric mean mass accretion rate Mg.m. as a value at 
which the corresponding DF of accretion rates reaches a 
maximum. All models in this section have a and 7 set 
to 10^^ and 1.4, respectively. Color symbols in Figure [5] 
present Mg.m. (in Mq yr~^) versus time-averaged stellar 
masses (M,) (in Mq) in the Class I stage of stellar evo- 
lution. The time-averaging is done over the duration of 
the Class I stage, which is distinct in each model. (In 
particular, model set 1 is plotted by red squares, model 
set 2 — by blue triangles-up, model set 3 — by green dia- 
monds, model set 4 — by black triangles-down, model set 
5 — by cyan circles, model set 6 — by red diamonds, model 
set 7 — by green triangles up, model set 8 — by blue dia- 
monds, and model set 9 — by black circles. Each filled 
symbol (of same color and shape) within a given set of 
models represents an individual object, which has formed 
from a cloud core of distinct mass, rotation rate, and size. 

In addition, each filled symbol is assigned vertical bars 
describing (for a given model) a typical range of mass 
accretion rates in the Class I stage. Defining these typical 
values turned out not trivial. This is because some of the 
models may show transient drops or rises in M, which are 
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Central object mass (M,.,) 

Fig. 5. — Geometric mean mass accretion rates versus time- 
averaged stellar masses in the Class I stage of stellar evolution. 
Symbols of specific shape and color represent different sets of mod- 
els as described in the text (see online version for the color figure). 
Vertical bars delineate typical variations in the mass accretion rate 
in every model. The lines are the least-squares best fits to all data 
(solid) and intermediate- and upper-mass objects (dashed). 



very short-lived and thus arc not statistically significant. 
Therefore, the vertical bars comprise only those accretion 
rates the DF for which is greater than 0.05. This means 
that we have excluded objects with extreme accretion 
rates, which have probability to be detected less than 
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one part in twenty. Considering the number of detected 
ob ject s along the li n e of c onstant central object mass (< 
20) in lEnoch et all (|2009D . we believe that our estimate 
represents conservative values. 

Figure [5] clearly demonstrates that the mass accretion 
rates in the embedded phase do not show any strong 
dependence on the stellar mass. The least-squares best 
fit to all data (solid line) yields the following relation 



0.25±0.05 



(13) 



which is considerably shallower than that for the TTSs 
and BDs. Accretion rates of the intermediate- and upper- 
mass objects (M, > 0.2 Mq) in Figure [5] seem to have a 
somewhat steeper dependence on (Af*), as indicated by 
the dashed line, yet the corresponding exponent 0.5 ±0.1 
is considerably smaller than that of the evolved TTSs 
and BDs of 0.5-3.0 Myr age. We again want to emphasize 
that our model does reproduce a steep dependence seen 
in TTSs and BDs for bo th purely self-gravitating disks 
(jVorobvov fc Basull2008D and self-grayitating disks with 
turbulent viscositv (jVorobvov fc Basij||2009bf ). 

A hint as to why this difference takes place can be 
seen in Figure [TJ The embedded phase is characterized 
by mass infall rates onto the disk Minfaii (dashed lines) 
that are often greater than the mass accretion rate onto 
the star (solid lines). This implies that accretion in the 
embedded phase may in part be controlled by the rate of 
gas deposition onto the disk, which in turn is expected 
to be only weakly dependent on the stellar mass. 

We demonstrate this in Figure[6l which shows the time- 
averaged infall rates (Minfaii) (filled circles) versus time- 
averaged stellar masses (M*) for most models in Table [TJ 
We have eliminated some of the models, which happen to 
have disk radii that exceed 600 AU — the radius at which 
we calculate the mass infall rate. The time-averaging is 
done over the duration of the Class I stage. For compari- 
son, we also plot the geometric mean mass accretion rates 

M (open circles) . Note that we use the same symbol type 
for all models in this figure. It is evident that (Minfaii) is 
systematically greater than Afg.m.. The least squares fit 
yields the following relation between the time-averaged 
infall rates and stellar masses 



(Afinfall) -4.0X 10-^(M,) 



0.20±0.03 



(14) 



the exponent of which is quite similar to that of the 
Afg.m. vs. (Af,) relation (dashed line, see also equa- 
tion [13]). A weak dependence of (Afinfaii) on the stel- 
lar mass is explained by a gradual increase in the enve- 
lope temperature along the sequence of increasing stellar 
masses. We also emphasize that our model infall rates 
lie wi thin the limits predicted by the standard model of 
IShui (jl977f) for star forming regions with gas temperature 
Tg = 10 - 25 K. 

7. CONCLUSIONS 

In this paper, we have studied mass accretion rates on 
to the star M and mass infall rates onto the disk Min{ai\ 
in the early embedded phase of stellar evolution, focusing 
on objects with stellar masses lying in the 0.08-2.0 Mq 
range. We have employed numerical hydrodynamic sim- 
ulations of the gravitational collapse of a large number of 
model cloud cores, which cover a wide range of masses, 
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Fig. 



- Geometric mean mass accretion rates onto the star 
(open squares) and time-averaged mass infall rates onto the 
disk (Afinfaii) (filled circles) versus time-averaged stellar masses in 
the Class I stage of stellar evolution. Most models from Table 1 
are shown in the figure. The dashed and solid lines are the least 
squares best fits to Mg.m. and (Minfaii), respectively. 



initial rotation rates, and sizes. The numerical integra- 
tion is started in the pre-stellar phase and followed into 
the late accretion phase when almost all of the initial 
cloud core has been accreted onto the star/disk system. 
We compare our model accretion rates with those re- 
cently inferred from observations of embedded YSOs in 
young star forming regions. We find the following. 

• Our numerical modeling yields a highly variable ac- 
cretion history in the phase that immediately fol- 
lows disk formation. Prolonged periods of low ac- 
cretion M < 10~^ Mq yr~^, when disks may be 
gravitationally unstable but are yet stable to frag- 
mentation, are interspersed with short-lived bursts 
of mass accretion with M > 10^^ Mq yr^-^, when 
disks forms fragments and drive them onto the cen- 
tral star via the exchange of angular momentum 
with spiral arms. The latter phenomenon is called 
the burs t mode of accretion and is investigated in 
detail in lVorobvov fc Basul (|2006[ j. 

• Model accretion rates M in the early embedded 
phase have a log-normal distribution, with its 
shape depending on the magnitude of turbulent vis- 
cosity (parameterized usin g the usual a-model of 
IShakura fc SunvaevI ()1973[ )) and disk temperature. 

• The spread in M (or the width of the accretion 
rate distribution function) is greater in models with 
lower a and smaller in models with higher a and 
higher disk temperature Td. An increase in either 
Td or a (or both) acts to stabiHze disk agains t 
gravitational instability (jVorobvov fc Basull2009aD . 
This suggests that gravitational instability may be 
a dominant cause of the accretion diversity in cir- 
cumstellar disks exposed to a continuing gas depo- 
sition from infalling envelopes. 

• Our models yield a large population of objects 
(37%-75%) with mass accretion rates M < 



10 
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10 ® Mq yr ^, which are smaller than those pre- 
dicted by the standard mode l of insid e-out col- 
lapse, a few X 10-6 Mq yr"! (lShulll977D . A simi- 
lar large fraction of embedded YSOs with sub-Shu 
accretion rates 55%) was recently reported by 
lEnoch et al.l (|200iD for Perseus, Serpens, and Ophi- 
uchus star forming regions. 

• Approximately l%-2% of objects in our mod- 
eling have super-Shu accretion rates M > 
10^^ Mq yr"^, which is roughly a factor of 2-3 
smaller than found by lEnoch et all (120091 ). This 
suggest that other (than disk fragmentation) mech- 
anisms may also contribute to bursts of mass ac- 
cretion in the embedded phase stellar evolution. 

• Purely viscous models, with disk self-gravity arti- 
ficially set to zero, significantly underestimate the 
number of objects with "standard" mass accretion 
rates lO"*' Mq yr"! <M< IQ'^ Mq yr-\ yield- 
ing essentially no objects with super-Shu rates. 
This demonstrates the importance of disk self- 
gravity in the early disk evolution. 

• Mean mass accretion rates in the early embed- 
ded phase show only a weak dependence on the 
time-averaged stellar masses of the form Afg.m. = 
10"6(M*)°-25±0 05, in sharp contrast to the steep 
(with exponent ~ 2.0) empirical relation inferred 
for evolved TTSs and BDs of 0.5-3.0 Myr age (see 



e.g. lMuzerolle et al.ll2003HNatta et al.ll2004) . This 
result is particularly interesting since our numeri- 
cal model can reproduce the steep relation in the 
late e volution phase for both p urely self-gravitating 
disks (jVorobvov fc Bas u 2008) a nd self-gravitatini 
disks with turbulent viscosity (jVorobyov &: Bas 
l2009bD . 



• The lack of strong dependence of Mg ,„ on (Af,) in 
the early embedded phase can be attributed to a 
continuing deposition of matter from the infalling 
envelope. The time-averaged rate of mass infall 
onto the disk (Afinfaii) is systematically greater 
(roughly by a factor of four) than the mean mass 
accretion rate onto the star Afg.„i., implying that 
A/g.m. is at list partly determined by (Minfaii). 
At the same time, (A/infau) exhibits only a week 
growth along the sequence of increasing stellar 
masses, which is caused by a moderate increase in 
the envelope temperature from ~ 10 K to ~ 25 K. 
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